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ABSTRACT 

As an initial investigation into the long-term evolution of protostcllar disks, we explore the conditions 
required to explain the large outbursts of disk accretion seen in some young stellar objects. We use one- 
dimensional time-dependent disk models with a phenomenological treatment of the magnetorotational 
instability (MRI) and gravitational torques to follow disk evolution over long timescales. Comparison 
with our previous two-dimensional disk model calculations (Zhu et al. 2009b, Z2009b) indicates that 
the neglect of radial effects and two-dimensional disk structure in the one-dimensional case makes only 
modest differences in the results; this allows us to use the simpler models to explore parameter space 
efficiently. We find that the mass infall rates typically estimated for low-mass protostars generally 
result in AU-scale disk accretion outbursts, as predicted by our previous analysis (Zhu et al. 2009a, 
Z2009a). We also confirm quasi-steady accretion behavior for high mass infall rates if the values of 
o-parameter for the magnetorotational instability is small, while at this high accretion rate convection 
from the thermal instability may lead to some variations. We further constrain the combinations of 
the o-parameter and the MRI critical temperature, which can reproduce observed outburst behavior. 
Our results suggest that dust sublimation may be connected with full activation of the MRI. This is 
consistent with the idea that small dust captures ions and electrons to suppress the MRI. In a later 
paper we will explore both long-term outburst and diskevolution with this model, allowing for infall 
from protostellar envelopes with differing angular momenta. 
Subject headings: accretion disks, stars: formation, stars: pre-main sequence 



1. INTRODUCTION 

In the standard model of low-mass star formation, 
a molecular cloud core collapses to a protostar over 
timescales of ~ 10 5 yr (e.g., Shu, Adams, & Lizano 1987), 
consistent with observations (Kenyon et al. 1990; Enoch 
et al. 2008). However, steady accretion of this mass onto 
central stars with a plausible mass-radius relationresults 
in accretion luminosities that are larger than those ob- 
served in low-mass protostars (Kenyon et al. 1990, 1994; 
Enoch et al. 2009). One solution to this "luminosity 
problem" is that most infalling matter first falls to the 
circumstellar disk and then is accreted to the star dur- 
ing short-lived outbursts; in this model protostars are 
usually observed in quiescence. The FU Orionis objects 
provide direct evidence for this type of behavior, with 
maximum accretion rates of 1O _4 M0 yr -1 over periods of 
decades to centuries (Hartmann & Kenyon 1996), which 
also directly suggests ~10 -2 M Q in the disk at ~ 1 AU. 

A number of theories have been proposed to explain 
FU Orionis outbursts, including thermal instability in 
the inner disk (Bell & Lin 1994), binary interactions 
(Bonnell & Bastien 1992), and gravitational clumping 
at several AU (Vorobyov & Basu 2005, 2006). In a re- 
cent paper (Z2009a), we explored the possibility that 
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outbursts might result because of a mismatch between 
the mass fluxes that can be transported by the mag- 
netorotational instability (MRI) in the inner disk, and 
the gravitational instability (GI) in the outer disk. Us- 
ing steady thin disk theory, we argued that outbursts 
are to be expected when disks are driven by mass addi- 
tion lower than 10 _4 MQyr _1 , as initially found by Ar- 
mitage, Livio, & Pringle (2001). We then developed a 
two-dimensional model of FU Orionis disks, which veri- 
fied that outbursts of accretion similar to those observed 
could be produced using reasonable parameters for MRI 
transport when thermal ionization dominates (Z2009 b). 

The computationally-intensive nature of two- 
dimensional (let alone three-dimensional) simulations 
of outbursts makes it difficult to conduct studies of the 
effects of differing parameters on disk evolution over 
significant timescales. We have therefore developed 
one-dimensional disk models to follow disk evolution. 
While such models have limitations, they can serve as 
a starting point to investigate the landscape of possible 
disk evolutionary pathways. The diversity of disk 
properties among stars of nearly the same age and mass 
(e.g., Hartmann 2009) is plausiblythe result of differing 
initial conditions, and models of the type we explore 
here can begin to address this possibility. 

In §2 we describe our 1-D, two-zone model, and com- 
pare its outburst properties with our 2-D models in §3. 
In §4 we present the results of a parameter study de- 
signed to show when outbursts occur and how the out- 
burst strength and frequency depend on the adopted pa- 
rameters. We present our discussion and conclusions in 
§5 and §6. The present work, in which we assume con- 
stant mass addition to the outer disk, serves as a starting 
point for our subsequent investigation of long-term disk 
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evolution with mass addition due to infall from a rotating 
protostellar envelope in a following paper. 

2. ONE-DIMENSIONAL TWO-ZONE MODELS 

In this paper we adopt a version of the "layered ac- 
cretion" disk model originally put forward by Gammie 
(1996). In this model, unless the disk is warm enough 
that thermal ionization is sufficient to couple the mag- 
netic field effectively to the neutral disk material, only 
an upper layer of the disk can sustain the MRI due to 
non-thermal ionization by cosmic and/or X-rays. A sig- 
nificant amount of work has been done to study the prop- 
erties of the active layer (e.g. Sano et al. 2000, Turner 
& Sano 2008, Bai & Goodman 2009), however due to 
the complex physics and chemistry involved, further the- 
oretical or even observation study is needed. Here, we 
assume that the mass column density which can be ion- 
ized is roughly constant with radius (see, e.g., Glassgold, 
Najita & Igea 2004). If the total disk surface density E 
in the "cold" regions is less than the limiting active layer 
column density E Q , the disk is assumed to be completely 
viscous with a given olm viscosity parameter due to MRI 
activity. On the other hand, if the disk surface density 
is higher than this limit, then only the surface layers are 
assumed to exhibit the MRI. 

Our one-dimensional, two-zone models (1D2Z) thus in 
general exhibit two layers: the surface layer and the cen- 
tral or "dead" zone 7 with surface density T, d . The dead 
zone is assumed to have little or no MRI activity, though 
it may exhibit transport due to GI (see below) . The tem- 
peratures T and T</ are averages which characterize the 
corresponding layers. 

The surface density evolves according to the mass 
conservation and angular momentum conservation equa- 
tions, 

2^ = ^, (1) 
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where Mj = —2TrT,iRvi is the radial mass flux in the disk, 
the stress Wr^a — KEiVidCl/ dR, and subscript i denotes 
either 'a' (active layer) or 'd' (dead zone). 

Assuming Keplerian rotation, equations (1) and (2) 
can be simplified to 

5 t E, = l dR (^d R (-R 2 W^)) , (3) 

where Or = d/dR. The constant infall rate M in is set 
as an inflow outer boundary condition with Mi„=10~ 4 , 
10~ 5 , and lO^Mgyr" 1 . 

The parameter E^ is the maximum non-thermally ion- 
ized surface density, assumed constant during each calcu- 
lation. If at a given timestep E a > E^, the excess mass 
of the active layer (Eq-E^) is added to the dead zone 
(Ed=Ed+E a -E / i). Conversely, if E a < Eyi and E^ ^ 0, 
part of the dead zone is assumed to be non-thermally 
ionized, in which case E a is set to be E^ and E^ de- 
creases to Ed — (E^ — E a ). Setting E^ = const, is a 

7 Here "dead" refers to magnetically dead without MRI but it 
could transport angular momentum due to the gravitational insta- 
bility. 



crude approximation; it is likely that E^ varies with ra- 
dius and depends on the local abundance of dust and flux 
of ionizing radiation. 

The temperatures are determined by the balance be- 
tween the heating and radiative cooling, 



CE,i9tTj — Qheat,i — Qcool,i 



(4) 



where the heat capacity is Cs^E^c 2 . ijTi. For the active 
layer, the cooling rate is determined by 
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where a is the Stefan-Boltzmann constant, T ext repre- 
sents the irradiation from the central star, and r a is the 
optical depth of the active layer. The final factor in equa- 
tion (5) is an approximate form which accommodates 
both optically thin and thick cooling. We assume 
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where L is the total luminosity of the star and f(R) 
accounts for the non-normal irradiation of the disk by 
the central star; here we set f(R) = const. = 0.1. The 
active layer optical depth is given by 



T a — ~^j(i K '{pcnT a ) . 



(7) 



where n is the Rosseland opacity derived from Z2009a at 
the active layer density and temperature, p a = Ti a /2H a , 
pd = (T, a + 'E d )/2H d , and H a and H d are the scale height 
of the active layer and the dead zone. 

The dead zone has a cooling rate similar to that of 
the active layer. If the active layer is optically thick, the 
incident radiation flux into the dead zone from the active 
layer is aT a . Thus the cooling rate is 



Qcooi,d — -^cr{T d 
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with 



Td = -T, d K(p d ,T d ) . 



(8) 
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However, if the active layer is optically thin, the incident 
flux becomes a(T a T a +T^ xt ), and 
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With increasing temperature the radiative cooling time, 
which constrains the numerical timestep, becomes small. 
For computational efficiency we then make the equilib- 
rium approximation 
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when the disk midplane temperature is > 5000K. 

The heating rate of the dead zone is just the viscous 
heating rate, while, in the active layer, the heating rate 
consists of its own viscous heating and the radiation from 
the underlying dead zone: 



Qheat,a — Qvisc ~t~ Qcool,d • 



(12) 
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In the case when the active layer is optically thin, this 
equation is modified to 



Qheat.a Qvisc T a Q cool,d • 

The viscous heating term is 

3 

Qvisc — 2 Wr^Q. , 



(13) 



(14) 



where Wr^ — (3/2)£jZ/jf2. To evaluate the viscosity z^, 
we have considered both MRI and GI transport. The net 
viscosity Vi is the sum of both, 



H 1T 



(15) 



where on — olq + olm and 



a Q = e~ Q . (16) 

The MRI viscosity is assumed to have a fixed value of 
olm whether the region in question is thermally or non- 
thermally ionized. We assume that above some critical 
temperature Tm the MRI is fully activated throughout 
the disk with viscosity parameter olm- The Toomre in- 
stability parameter Q is evaluated using the disk cen- 
tral (midplane) temperature T^, the total surface density 
(Ea+Sd), and assuming Kcplerian rotation. The form of 
oiq is motivated by a desire to make gravitational torques 
significant only when Q < 1.4, as indicated by global 
three-dimensional simulations (e.g., Bolcy ct al 2006). 

There are uncertainties in adopting the above approach 
to transport. GIs involve large scale density waves that 
cannot be captured with a local viscous treatment, al- 
though local treatments are adequate under some cir- 
cumstances (Lodato & Rice 2004; Cossins et al. 2009). 
As discussed in Z2009a, the essential properties of this 
treatment are the assumptions that disks with GI have 
Q-values of order unity, and that the GI produce local 
dissipation of the accretion energy Under these assump- 
tions, the precise form of olq will not affect the disk's 
evolution, as long as a is a steeply declining function of 
Q near ~1.5. To make this point clear, we ran the same 
simulation for a test case but with the olq prescription 
of Lin & Pringle (1987,1990). As expected, the different 
forms of «q have no effect on the disk outbursts. 

Similarly, whether the MRI can be fully activated 
in a non-thermally-ionized layer depends in part upon 
whether small dust grains have been sufficiently depleted 
(e.g., Sano et al. 2000). This is a complicated problem 
with substantial observational and theoretical uncertain- 
ties; we therefore adopt the simplest possible approach. 
It turns out that the value of E a is unimportant for un- 
derstanding large outbursts (Z2009a), as long as the GI 
in the dead zone transports more mass than the active 
layer; but E a does have important effects on the long- 
term disk evolution at low accretion rates, as discussed 
in a following paper. 

It is now clear that the magnetic fields that give rise to 
ajf diffuse radially (Lesur & Longaretti 2008, Guan & 
Gammic 2009, Fromang & Stone 2009) and take time to 
build up and decay (e.g. Hirose et al. 2009). To account 
for these effects we introduce an evolution equation for 
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where olm,o is the equilibrium value for olm- The first 
term permits oim to relax up, or down, as Tm is crossed. 
The second term corresponds to radial diffusion of the 
magnetic field. For numerical reasons we set the dimcn- 
sionless radial diffusion coefficient to 0.5 (the radial dif- 
fusion coefficient is actually a function of distance from 
the midplane). 

3. OUTBURST BEHAVIOR 

In the protostellar phase, the disk is unlikely to trans- 
port mass steadily from ~100 AU all the way to the 
star at an accretion rate matching the mass infall rate 
10~ 6 — 1O~ 4 M yr _1 from the envelope to the outer disk. 
This mismatch leads to outbursts which are qualitatively 
similar to that found by Armitagc et al. (2001), Book & 
Hartmann (2005), and in our 2-D hydrodynamic simula- 
tions (Z2009b). In summary, before the outburst, mass 
added to the outer disk moves inwards due to GI, but 
piles up in the inner disk as GI becomes less effective at 
smaller radii. Eventually, the large E and energy dis- 
sipation leads to enough thermal ionization to trigger 
the MRI at several AU. The MRI front quickly moves 
in across the inner disk and the inner disk accretes at a 
higher mass accretion rate, resembling FU Orionis-type 
outbursts. This high mass accretion rate during the out- 
burst also makes the inner disk thermally unstable. After 
the inner disk has been drained by the outburst and be- 
comes too cold to sustain the MRI, the disk returns to 
the low state. With the mass continuously accreted from 
the outer radii (or from an infalling envelope), the disk 
evolves to conditions leading to another outburst. This 
MRI triggered by GI outburst can also be understood 
as the classical thermal instability, but the S-curve is 
formed primarily by the variation of a near Tm, rather 
than variations in opacity and assumed variation in a 
near hydrogen ionization (e.g. Bell & Lin (1994)). 

To test how 1D2Z models simulate the outbursts com- 
pared with 2-D simulations, we set up a test case with 
all the parameters adopted from our previous 2-D sim- 
ulations (Z2009b). In both 1-D and 2-D simulations, 
we have used an updated opacity from Z2009a and 
Tm=1500 K . Because Z2009b do not consider irradia- 
tion, the irradiation factor f in equation (6) for the 1D2Z 
simulation is set to be 0. The inner radius in both cases 
is set to 0.2 AU. 

Figure 1 shows the mass accretion rate as a function 
of time for both 1-D and 2-D simulations. As shown, the 
1-D simulations closely resemble 2-D simulations at the 
equilibrium states, such as the state before the outburst 
is triggered and the state during the outburst. For some 
rapid, or small scale, disk variations, such as the MRI 
front propagation and the convective eddies in the hot 
inner disks, the 2-D simulations exhibit more complex 
behavior than the 1-D simulations, so the outbursts dif- 
fer in detail. In particular the 1-D simulations show an 
initial high M peak at the beginning of the outburst that 
is not seen in 2-D. The 1-D simulations also show "drop 
outs" in accretion that do not occur in 2-D. 

In detail, starting from the MRI activation at ^2 AU, 
the MRI active region moves inwards. During this pro- 
cess, mass piles up at the inner boundary of the active re- 
gion, because the disk is MRI active beyond this bound- 
ary and has a higher mass accretion rate than the MRI 
inactive disk at smaller radius. In 1-D simulations, which 
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do not capture the effects of radial pressure gradients, 
mass piles up at this boundary. When this mass eventu- 
ally accretes on to the central star there is a sharp peak 
in M. 

The M drop-outs during outbursts in the 1-D simu- 
lations are related to the thermal instability associated 
with hydrogen ionization. The drop-outs occur when the 
inner disk returns to the TI low state in 1-D; in 2-D ra- 
dial and vertical convection smooths out variations in the 
accretion rate and allows the inner disk to remain in the 
high state. 

Despite the initial M peak and the drop-outs in the 
1-D model, the outburst timescale, M, and consequently 
the total mass accreted during one outburst are similar in 
1-D and 2-D simulations (Fig. 1). This similarity is due 

to the fact that the outburst timescale and M are just 
determined by the radius where the MRI is triggered and 
by cum! the total mass accreted during one outburst is 
the mass difference before and after the outbursts, which 
are both equilibrium states. 

4. PARAMETER STUDY 
Having tested that the 1-D models reproduce the gen- 
eral properties of the outbursts (maximum M, duration 
time, etc.), we next turn to a parameter study with 1- 
D models to test the predictions of Z2009a using steady 
state models. 

The radial range considered is from ^0.2 AU to 30 
AU. We adopt a constant inflow boundary condition at 
30 AU. The inner boundary 0.2-0.3 AU is chosen to avoid 
instabilities which occur at the smallest radius where 
there is a transition between thermal and non-thermal 
activation of the MRI, here called Ri. 

We find that, at least for our phenomenological model, 
Ri is unstable. Wiinsch et al. (2006) has also found 
even with a 2-D radiative transfer hydrodynamic layered 
model, this instability still occurs. Ri oscillates around 
a mean radius, with thermal fronts washing inward and 
outward. This is similar to thermal instability but due to 
a variations between the dead zone and the inner MRI 
active region (Wiinsch et al. 2006). However, we should 
be cautious in this instability, especially since dust may 
be sublimated before MRI activation. Then the inner 
disk may become optically thin due to the low opacity of 
the gas, and direct irradiation from the central star would 
ionize the dust wall. A proper treatment of Ri needs a 
3-D MHD simulation with irradiation, dust physics, and 
ionization physics. Because this treatment is impracti- 
cal, we set the inner boundary of our model just outside 
Ri. Ri depends on a variety of parameters (central star 
luminosity, the irradiation angle to the disk's surface, ac- 
tive layer surface density, the MRI trigger temperature 
and its viscosity parameter). Generally, Ri increases as 
the heating rate increases. In the end we set the inner 
boundary to 0.2AU if a M =0.01 and 0.3 AU if a M =0.1. 

The infall is treated as a constant mass inflow at the 
disk outer edge at 30 AU. This assumption may not be 
applicable to study the disk's long-term evolution during 
the entire infall stage (10 5 years) since the infall centrifu- 
gal radius increases with time; this is explored exten- 
sively in Paper II (Zhu et al. 2010b). Over the timescale 
of a single outburst, as studied in this paper, this as- 
sumption is a good approximation as long as the infall 



centrifugal radius is larger than the MRI stable radius 
(R M 1-10 AU in Figure 8). 

Because protostars are thought to form over ~ few 
xlO 5 yr, the average infall rate should be 10~ 6 - 
lO^Moyr" 1 (Stahler 1988; Hartmann et al. 1997); nu- 
merical simulations suggest the infall rate could be up 
to 10~ 4 M Q yr^ 1 at the earliest stages (Bate et al. 2003). 
Thus we study the disk evolution with infall rates varying 
from 10~ 4 M Q yr" 1 to 10~ 6 M Q yr" 1 . 

We neglect any change in the central star mass over 
the outburst timescale, which is a reasonable approxi- 
mation. The effect of changing central mass over longer 
evolutionary timescales is included in Paper II. 

What is cum? Observations of dwarf novae, X-ray bi- 
naries, and FU Ori suggest a ~ 0.1 (King et al. 2007; 
Zhu et al. 2007). T Tauri disk observations suggest 
a ~ 0.01 (Hartmann et al. 1998). MHD simulations sug- 
gest a > 0.01, with the precise value depending on the 
resolution, treatment of small-scale dissipation, stratifi- 
cation, and treatment of radiation transport (e.g. Fro- 
mang & Papaloizou 2007, Guan et al. 2009, Davis et al. 
2009, Shi et al. 2009; Hartmann et al. 1998). Because 
the situation is not yet resolved, we will consider cases 
with olm — 0.01 and om =0.1. 

The MRI activation temperature Tm is not known pre- 
cisely (it depends on the location of alkali metals, their 
ionization rate, the abundance of small grains, and the 
threshold ionization fraction for MRI turbulence), so we 
consider cases with Tm=1400 K and 1800 K. These val- 
ues are chosen, consistent with our opacity prescription, 
to represent cases with and without dust. 

4.1. Dependence on M in 

Figures 2 - 4show the outbursts with different infall 
rates varying from 10 _4 M© yr -1 to 1O _6 M0 yr -1 . For a 
given set of disk parameters («m, Tm), the disk accretes 
quasi-steadily if the infall rate is high enough (small vari- 
ations could still appear due to the classical thermal in- 
stability, as discussed in §5.1). If au — 0.01, the disk 
accretes steadily if the infall rate is 10 _4 M Q yr _1 , while 
outbursts appear with smaller infall rates. With a bigger 
olm = 0.1, the disk accretes nonsteadily/in outburst for 
all the cases with infall rates 10" 4 M Q yr- 1 -10 _6 M Q yr" 1 . 
This is consistent with our predictions from the steady 
state models (Z2009a). 

.The outbursts become shorter with smaller infall rates 
Mi n . For example, with om=0.1 and Tm=1400 K, the 
outbursts last 800, 500, and 400 years for M in = 10~ 4 , 
10~ 5 , and 10~ 6 M Q yr _1 .This can be explained by the 
fact that the MRI is triggered at larger radii with larger 
M in and thus the outburst (viscous) timescale is longer 
(§5.1). The simple analytical calculation in the Appendix 
shows that the outburst timescale is 

(18) 

for T M = 1400 K, where Rq is the Gl-induccd MRI acti- 
vation radius which will be discussed in §5. This agrees 
with the numerical simulations reasonably well. The 
mass accretion rates during the outbursts are similar. 
The infall rate also determines the outburst fre- 
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quency. The time between two outbursts issignificantly 
shorter with M in = lO _4 M yr _1 than withM in = 
10 _5 M Q yr _1 , because the mass accreted to the star 
during all the outbursts should be equal to the mass 
from infall integrated over the same period of time. 
Since the outburst timescale is insensitive to the infall 

• 1/9 

rate (~ Mj„ ), we assume each outburst transports 
~ O.O3M , which is suggested by the observation of FU 
Ori, and thus the timescale between two outbursts is 0.03 
Mq I Mi n ■ Therefore higher infall rates lead to more fre- 
quent outbursts. 

4.2. Dependence on cum and Tm 

The effect of qm on the outburst can be seen by com- 
paring the upper and lower panels ofFigures 2-4. For 
a given Mi n and Tm, with a higher ajf the outburst 
is shorter and stronger. This is because the outburst 
timescale is close to the viscous timescale, which is in- 
versely proportional to a M while the mass accretion rate 
is proportional to aM- This can be understood using 
Equation (18) and 

M = 5x lO-4^M yr- 1 , (19) 

as shown in the Appendix for the Tm=1400 K case. 

On the other hand, comparing the left and right panels 
of Figures 2 - 4, we find that the disks with higher Tm 
have shorter but stronger outbursts. The outburst is 
shorter because the MRI is triggered at a smaller radius 
if Tm is higher, resulting in a shorter viscous timescale. 
Figure 5 shows the disk's condition just before the MRI 
is triggered in the case Tm=1400 K (the solid curve) and 
7m =1800 K (the dotted curve). The MRI is triggered 
at 3 AU with T M =1400 K and 1.5 AU with T M =1800 
K. The outburst is stronger with higher Tm because the 
surface density at the MRI trigger radius with Tm= 1800 
K (1 AU) is much higher than the surface density at the 
MRI trigger radius with T M =1400 K (3 AU) (upper left 
panel in Fig. 5). Thus, with higher Tm, the smaller but 
more massive inner disk leads to shorter but stronger 
outbursts. 

Our simulations also indicatethat the outbursts with 
higher Tm accrete less mass than those with lower Tm 
no matter what the infall rate is (Table 1), because the 
outbursting inner disk extends to smaller radii, and thus 
contains less mass, for a higher Tm- 

4.3. Dependence on M* and Ha 

Protostars have a variety of masses, most of which 
are less than 1 M Q . We also calculated results for 
central star masses ( 0.3, and 0.1 M Q ) with the same 
T M (1400 K) and a M =0.1, and the infall rate M in = 
10~ 5 M Q yr _1 .The mass accretion rates with time for 
these cases are shown in Figure 7. The outbursts have 
similar mass accretion rates, but the outburst is slightly 
shorter with a less massive central star, as suggested by 
equation (18). In addition the mass accreted during an 
outburst is not significantly affected by its central star 
mass (Table 1). The outbursts are broadly similar, even 
though the central star masses differ by a factor of 10. 

Another uncertainty in a layered disk model is ZU, 
which depends on the flux of ionizing radiation and 



the abundance and size distribution of the dust. In 
protostars with mass infall rates from lO _6 M0yr _1 to 
10~ 4 M Q yr _1 , most of the infall mass is transported to 
the inner disk by the GI in the dead zone, and mass trans- 
port through the active layer is negligible. Thus XU has 
little effect on the unsteady accretion and the outburst 
(Fig. 6), but it does affect the mass accretion rate in 
the low state. The mass accretion rate in a layered disk 
is determined by the active layer mass accretion rate at 
Ri. At R < Ri the disk is MRI active due to thermal 
ionization. Thus the disk mass accretion rate should be 
proportional to Hacxm, which is shown in Figure 6. 

5. DISCUSSION 
5.1. Unsteady accretion region 

The above parameter study can be simply summa- 
rized in the M — R plane as shown in Figure 8, which 
builds upon the steady state vertical structure calcula- 
tions (Z2009a; an analytical analysis is given in the Ap- 
pendix). The shaded regions in Figure 8 are the radii 
at which disk accretion is expected to be unsteady if the 
mass infall rate is constant. 

The solid curve farthest to the lower right in the Figure 
(labeled Rq) is the radius where the central temperature 
of a pure Gl-driven disk (Q=l disk) would reach Tm (e.g. 
1400 K). In other words, at a given infall rate Mj„, pure 
GI disks can accrete steadily beyond Rq by the GI, but 
within Rq the MRI will be activated. Up and to the 
left of Rq in Figure 8 another solid curve, labeled Rm, 
denotes the radii at which a pure MRI disk of the given 
aM would have a central temperature of Tm- If the disk 
is MRI active, it can accrete steadily purely by the MRI 
within Rm, but beyond Rm the MRI will be turned off. 
When Rm and Rq cross a smooth transition between 
the GI and MRI exists and steady accretion is possible. 
From the left panel of Figure 8, we see thatthe disk can 
accrete steadily with a=0.01 and Mj„ = 10 _4 M Q yr~ 1 . 
This is observed in the time evolution discussed in §4.1 
and Figure 2 (see also Armitage et al. (2001)). 

We predict that in the shaded regions matter will pile 
up through the action of GI, trigger the MRI, and pro- 
duce an outburst. The dotted curve shows Rq and Rm 
if Tm=1800 K (at which temperature all dust has sub- 
limated). Rq and Rm at 1800 K arc smaller than they 
are at 1400 K because of the temperature plateau around 
the dust sublimation temperature (Z2009a; also can be 
seen at R- 2-10 AU in Figs. 5). Thus if T M is higher, 
outbursts are expected to be shorter because the out- 
burst drains the small inner disk (R < Rq) on a shorter 
viscous timescale (as demonstrated in Figures 2-4). The 
effect of aM can also be seen by comparing the left and 
right panels of Figure 8. 

The classical thermal instability will also be triggered 
at the infall phase as shown by the upper left shaded 
band in Figure 8. Even the 'steady'accretion case dis- 
cussed above with cn=0.01 and M in = 10 _4 M Q yr _1 is 
subject to thermal instability at R < 1AU. The two lines 
shown in Figure 8 correspond to the two limiting values 
of the "S curve" (e.g., Faulkner, Lin, & Papaloizou 1983) 
at which transitionsto the high (rapid accretion) state 
and the low (slow accretion) state occur. The instabil- 
ity depends on the disk's vertical structure (different "S 
curves") which behavesdifferently in 2-D than in 1-D, as 
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discussed in §3. We expect nature to behave somewhat 
more like the 2-D than thel-D case, so for the steady ac- 
cretion model discussed above convection will add small 
variations in M andwe refer to these steady accretion 
cases as quasi-steady. Generally, the thermally unstable 
region is distinctive at M > 10 -5 M© yr" -1 , which sug- 
gests the TI may be common for protostellar disks. 

The solid dots in Figure 8 are the MRI trigger radii 
from our 1-D simulations when Tm — 1400K. Although 
the trigger radii do not fall precisely on Rq, most of 
them are in the shaded region, indicating that this M — R 
plane has predicted non-steady accretion, with potential 
outbursts, to occur for infall rates < lO _5 M yr _1 for 
a M = 0.01 and < lO^Moyr -1 for a M = 0.1. For M = 
lO^Mgyr" 1 and a M = 6.01, the MRI-GI instability 
won't occur based on this M — R plane, which agrees 
with the 1-D simulations. Generally, our M — R plane 
results provide a good guide to the parameters for which 
unsteady accretion occurs. 

The outbursts and TI inquasi-steady accretion could 
provide a much hotter thermal history for the protostel- 
lar disks, which may have imprints on meteorites and the 
chemicals in the protostellar disks. 

5.2. Steady vs Non-steady accretion 

The non-steady disk accretion in our model is the result 
of the inability of the inner disk to transport mass inward 
at the same rate as mass is fed in from the outer disk by 
infall, 1O~ 5 M yr -1 . This mass pileup eventually leads 
to MRI activation and outbursts. Even after infall ends, 
the dead zone evolves with time if the mass addition from 
the outer disk due to GI and viscous stresses exceeds 
what can be carried inward by the active layer. 

Terquem (2008) was able to construct steady-state disk 
accretion solutions with active layers and "dead" zones. 
The difference is due to the assumptions of both finite, 
non-GI viscosity in the dead zone and to assuming a 
very much lower mass accretion rate, 10 _8 M Q yr _1 . At 
such low accretion rates, even very low non-GI dead zone 
viscosities and comparable active layer properties suffice 
to transport mass at these rates through the inner disk, 
which is not possible in our disks driven on the outside at 
high mass infall rates (except for small centrifugal radii, 
in which case the disk becomes fully MRI-active) . 

5.3. Constraints from observations of FU Orionis 
objects 

FU Orionis objects are outbursting systems withmax- 
imum disk accretion rates M ~ lO _4 M0yr _1 and a de- 
cay time of decades to hundreds of years (Hartmann & 
Kenyon 1996). These properties constrain our parameter 
space. 

The outbursts produced by 1D2Z models are sensi- 
tive to aM and Tm- Figures 2-4 show that qm=0.1, 
Tm=1800 K leads to outbursts that are too strong 
(M=10" 3 M Q yr" 1 ), while a M =0.01, T M =1400 K leads 
to outbursts that arc too long (3000 yr) 8 . If T M =1400 
K, the outburst start at ~3 AU (Figure 5), so aM needs 

8 The observed decay time is = (din F/dt)^ 1 (F is the flux; 
thus Td is an e-folding time). In principle the outburst could have 
noncxponcntial time dependence and the duration of the outburst 
could differ from r<j. In our models, however, the luminosity ex- 



to be large to produce the correct decay timescale. If 
Tm = 1800K the outburst is triggered at ~ 1AU, so 
«m ~ 0.01 is required to maintain the correct outburst 
timescale. 

One uncertainty here is that qm in a non-thermally 
ionized active layer may not be the same as om in a 
thermally ionized, dust free, fully conducting medium. 
Thus we mostly constrain qm in the latter (outburst 
stage) from the decay timescale; aM in the active layer 
is not well-constrained by outbursts. 

Another uncertainty is related to our assumption of 
constant Ha- If is instead a function of radius the 
disk's long term evolution will change. At early stages, 
however, if the GI in the outer disk dominates the disk's 
accretion, variation of XU will have little effect on the 
outburst mechanism discussed here. 

The central star mass and infall rate do not signifi- 
cantly change the outbursts. This agrees with the ob- 
servation that FU Orionis outbursts occur for protostars 
with different infall rates (Quanz et al. 2007 and Zhu ct 
al 2008 found that FU Orionis objects could be either 
Class or Class I objects.). 

Zhu ct al. (2007) argued that the decay timescale for 
FU Ori implies a ~ 0.02-0.2. This conclusion can be 
tested with our 1-D simulations. The last column in 
Table 1 shows the viscous timescale calculated using Rm 
and T M - 



where 

"=4> (2i) 

c s and n are calculated with T M and R M . 

Comparing with the outbursts' duration in column 6, 
we see the viscous time is close to the outburst time 
for T M =1400 K. If T M = 1800 K, the viscous time is 2- 
3 times longer than the outburst time, which may be 
due to the disk's temperature during the outburst being 
higher than T M if the MRI is triggered at smaller radius 
as Tm=1800 K case. However, since constraining a by 
using the viscous timescale is anorder of magnitude esti- 
mate, our simulations are consistent with the Zhu et al 
(2007) estimate. 

6. CONCLUSIONS 

In this paper, we have evolved a one-dimensional lay- 
ered disk model including both MRI and GI to study the 
unsteady disk accretion of protostars. The 1-D models 
reproduce the general properties of 2-D (axisymmetric) 
outbursts reasonably well, such as the outbursting mass 
accretion rate, durationand the accreted mass during one 
outburst. Because the 1-D model is faster, itenables us 
to study outbursts in an extended parameter space. 

Our results confirm that the disk can accrete steadily 
with high infall rates (M m ~10- 4 M Q yr" 1 if a M =0.01; 
Armitage et al. (2001)). This steady accretion may still 
have short timescale variations, however, due to the ther- 
mal instability in the inner disk, as suggested by our ear- 
lier, 2-D simulations. 

ceeds the preoutburst luminosity for a time only slightly longer 
than Trf. 
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We also confirm the prediction by Z2009a that pro- 
tostars are likely to accrete unsteadily /in outbursts for 
infall rates < lO^Moyr" 1 with a M = 0.01 and < 
10 _4 M Q yr _1 for am = 0.1. Outbursts are triggered 
at r ~ 1 — 10 AU for protostellar infall rates ~ 10~ 5 — 
10 M0JT -1 . The outbursts are stronger and shorter 
with larger c*m or Taj. The total mass accreted dur- 
ing one outburst mainly depends on Taj. While the 
outbursts are slightly shorter for more massive central 
stars, the outburst M is nearly independent of central 
star mass. The active layer surface density only affects 
the mass accretion rate inthe low state; it has little effect 
on the outburst. 

By comparing with the mass accretion rate and dura- 
tion ofobserved FU Orionis events, we can constrain a 
combination of «m and Tm- If «m is low (0.01), Tm 



needs to be high (1800 K, higher than the dust sublima- 
tion temperature); if q« is high (0.1) then Tm needs to 
be low (<1400 K). 

Our results show that 1-D, two zone models can cap- 
ture the basic features of disk evolution, given our as- 
sumptions about the action of the MRI and GI. In a later 
paper we will address disk evolution over a much longer 
timescale, explicitly taking into account mass infall from 
a rotating protostellar cloud. 
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APPENDIX 



By assuming a marginally gravitationally stable (Q=1.5) disk, the disk's structure is determined with a given M, 
and thus the radius where the outburst is triggered (Rq in Fig. 8 and Z2009a) can be derived. Unlike Z2009a where 
the detailed vertical structure is calculated numerically to give Rq, here we give simple analytical results for Rq by 
assuming the disk is vertically isothermal with constant opacity at a given radius. 

First, if k=CT q P /3 , the relationship between E and the central (midplane) temperature T c is given by 



^ = - 8 T* ff T=^T* ff X K . 
Using the form for k and using i o c =E/2iJ, where H is the disk scale height, 

T =3 1 /(4-a-/3/2)2(-4-/3)/(4-a-/3/2)ji4/(4-a-/3/2) s (i +i g)/(4-a-/3/2) c l/(4-a-/3/2) 



(1) 



xJ yj/(4-a-/V2) 



k s /3/(8-2a-/3) 

[im H J 



or equivalently 



s = 3 -l/(l +/ 8) 2 (4+/3)/(l+/3) T -4/(l+/3) T (4-a-/3/2)/(l+/3) (7 -l/(l+«Q-/3/(l+/3) 

e ff c 



-/3/(2+2/3) 



(2) 



(3) 



where f2 is the angular velocity at R, k is the Boltzmann constant and m# is the unit molecular mass. 

Then with Q=c,,f2/7rG£, and inserting equation (3) into Q to derive the relationship between R and T c at a given 
Q and M 

2/(9+6/3) 



jR = 3 4/(9+6/3) 2 (-14-2/3)/(9+6/3) 7r (-4-2/3)/(9+6/3) 



Mi, 



j.(2a+2/3-7)/(9+6,3)£2/(9+6/3) 



x £,1/(9+6/3) M (3+2/3)/(9+6/3) 



\ (l+2/3)/(9+6/3) 



Q 



(-2-2/3)/(9+6/3) 



(4) 



The dust opacity fitting from Z2009a suggests, at T< 1400 K, a=0.738, /3=0 and C=0.053. If we plot the relationship 
between R and M by given T c =T Af =1400 K and Q=1.5, we find 

lO- 4 M yr-i) \M Q , 



Rq = 11 



AU. 



(5) 



which corresponds well with the Rq calculated in Figure 8. 

bince RqocM 2 / 9 for T M =1400 K case, the outburst timescale is roughly 



R Q 



960 



Mi, 



0.1 



a M \ 10- 4 M Q yr- 



1/9 



M 



2/3 



yr, 



(6) 



where v is calculated for T c =Tm=1400 K. During outburst, however, T c could be higher than the MRI trigger 
temperature Tm, especially in the inner part of the disk. 
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If we further assume the outburst M ~ i/H(Rq) (this is the steady accretion disk solution which may not be true in 
the time-dependent case), we find 

M~3wX~ 3aM J° Cst . (7) 
GQ 

Here we explicitly distinguish between the sound speed c so during the outburst and the sound speed c st before MRI 
activation. For an order of magnitude estimate we assume c so ~ c st ~ c s (1400K) and so 

M^SxlO-^Moyr- 1 (8) 

which agrees with the numerical simulations for T M =1400 K reasonably well considering we are using a steady state 
assumption. 

We have only applied the above equations for theTM=1400 K case, since the detailed vertical structure is important 
if Tm=1800 K where the midplane is dust-free while the surface has dust. Also the outburst is triggered at a smaller 
radius for Tm=1800 K, and thus the disk temperature during the outburst can be much higher than the MRI trigger 
temperature. 
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Fig. 1. — The mass accretion rate with time for both 1-D (dotted curve) and 2-D (solid curve) simulations. 
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Fig. 2. — The mass accretion rate with time for different cxm and Tm for the mass infall rate of 1O -4 M0 yr _1 . 
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Fig. 3. — The mass accretion rate with time for different olm and Tjj for the mass infall rate of 10 _5 Mq yr~ 1 . Compared with the solid 
curves which are from the simulation with OQ=exp(-Q 4 ), the dotted curve in the lower left panel shows the outburst obtained using the 
qq prescription of Armitage et al. (2001). 
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Fig. 5. — The disk's radial structure at the stage just before the MRI is triggered in cases where Tjf =1800 K (dotted curve) and Tm=1400 
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14 



Zhu, Hartmann, & Gammic 



-6 



o 



8 








500 1000 1500 2000 

t(yr) 



Fig. 6. — The disk mass accretion rate with time for different active layer surface density (50 g cm 2 for the dotted curve and fO g cm 2 
for the solid curve). «^ =0.1 and Tm=1400 K. 
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Fig. 7. — The mass accretion rate with time for different central star masses: 1 Mq (solid curve), 0.3 Mg (dotted curve), 0.1 Mq (dashed 
curve). The infall rate is 1O~ 5 M yr" 1 , T M =1400 K, and a=0.1. 




Fig. 8. — Unstable regions in the R — M plane for a IMq central star. The shaded region in the lower right shows the MRI-GI instability 
with the MRI trigger temperature of 1400 K. The dotted curves show Rm and Rq (the boundaries of the MRI-GI instability shaded region; 
see text for definition) for an MRI trigger temperature of 1800 K. The shaded region in the upper left shows the region subject to classical 
thermal instability. The solid dots represent the radii where the MRI is triggered in 1-D simulations for Tm=1400 K. 
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TABLE 1 

1D2Z MODELS 



M* 


infall rate 


a 


Tm 


outburst M a 


duration 


accreted mass b 


Rm c 


viscous timc d 


M 


Moyr" 1 




K 


Mq yr" 1 


yr 


Mq 


AU 


yr 


1 


10- 4 


0.1 


1400 


2xl0~ 4 


700 


0.7 


6 


428 


1 


10- 4 


0.1 


1800 


10~ 3 


50 


0.02 


1.8 


182 


1 


10~ 5 


0.1 


1400 


2xl0~ 4 


400 


0.057 


3.5 


327 


1 


10~ 5 


0.1 


1800 


10~ 3 


40 


0.024 


1.2 


150 


1 


10~ 5 


0.01 


1400 


5xl0~ 5 


4000 


0.045 


3.2 


3127 


1 


io~ 5 


0.01 


1800 


8xl0~ 5 


400 


0.01 


1.2 


1490 


1 


10~ 6 


0.1 


1400 


2xl0~ 4 


350 


0.08 


1.9 


240 


1 


10~ 6 


0.1 


1800 


10~ 3 


50 


0.02 


0.85 


125 


1 


10~ 6 


0.01 


1400 


5xl0~ 5 


3000 


0.04 


2.1 


2533 


1 


10~ 6 


0.01 


1800 


8xl0~ 5 


800 


0.015 


0.87 


1268 


0.3 


10~ 5 


0.1 


1400 


2xl0~ 4 


300 


0.04 


2.3 


145 


0.1 


10~ 5 


0.1 


1400 


2xl0~ 4 


250 


0.025 


1.43 


66 



a the mass accretion rate at the half-time of the outburst 
b the mass accreted during one outburst 
c the MRI trigger radius 

d Thc viscous timcscalc is calculated by using Rm and Tm- 



